/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package javalens.astronomy;

import javalens.util.Position;
import javalens.util.ComplexNumber;
import javalens.util.CalculationException;
import javalens.util.Util;
import java.util.ArrayList;
import java.util.List;
import javalens.globalconfig.config;
import javalens.maths.MathsEngineFactory;

/** Basic lensing calculations. WORK IN PROGRESS
 *
 * @author vermaak
 */
public class Calculations {

    public static double calcSingleLensAmplification(Event event, double day) {
        // first obtain position in source plane. 
        Position sourcePos = calcSourcePlanePosition(event, day);
        double x = sourcePos.getX();
        double y = sourcePos.getY();
        double u = Math.sqrt(x * x + y * y);
        double amp = (u * u + 2.0) / u / Math.sqrt(u * u + 4.0);
        return blend(event.getF(), amp);
    }

    public static Position calcSourcePlanePosition(Event event, double day) {
        final double lambda = 4.726;
        final double tp = 1099.0;
        final double beta = -0.08;
        final double e = 0.0167;

        // psi in radians
        double u, T, funnyphi, tptime, p0, sqgl, aries;
        double poft, doft, x1, x1tmax, x2, x2tmax, start;
        double radang;
        double cosradangminusninety, sinradangminusninety;

        aries = 77. * Math.PI / 180;
        funnyphi = aries + lambda + Math.PI;
        T = 365.26;

        double tm = event.getTm();
        double rho = event.getRho();
        double te = event.getTe();
        double psi = event.getPsi();
        double b = event.getB();
        double ang = event.getAng();

        // calculate x1tmax and x2tmax
        tptime = (2 * Math.PI * (tm - tp) / T);
        sqgl = tptime + 2 * e * Math.sin(tptime) - funnyphi;

        x1tmax = rho * (-Math.sin(beta) * Math.cos(funnyphi) *
                (Math.cos(sqgl) - e) - Math.sin(beta) * Math.sin(funnyphi) * Math.sqrt(1. - e * e) * Math.sin(sqgl));
        x2tmax = rho * (-Math.sin(funnyphi) * (Math.cos(sqgl) - e) + Math.cos(funnyphi) *
                (1. - e * e) * Math.sin(sqgl));

        // Now for today's point
        tptime = (2 * Math.PI * (day - tp) / T);
        p0 = (day - tm) / te;
        sqgl = tptime + 2 * e * Math.sin(tptime) - funnyphi;

        x1 = rho * (-Math.sin(beta) * Math.cos(funnyphi) *
                (Math.cos(sqgl) - e) - Math.sin(beta) * Math.sin(funnyphi) * Math.sqrt(1. - e * e) * Math.sin(sqgl));
        x2 = rho * (-Math.sin(funnyphi) * (Math.cos(sqgl) - e) + Math.cos(funnyphi) *
                (1. - e * e) * Math.sin(sqgl));

        poft = p0 + Math.cos(psi) * (x1 - x1tmax) + Math.sin(psi) * (x2 - x2tmax);
        doft = b - Math.sin(psi) * (x1 - x1tmax) + Math.cos(psi) * (x2 - x2tmax);

        // with prlx done, we convert poft and doft to the real values in the plane
        radang = ang * Math.PI / 180.;
        cosradangminusninety = Math.cos(radang - Math.PI / 2.0);
        sinradangminusninety = Math.sin(radang - Math.PI / 2.0);

        double sx = poft * cosradangminusninety - doft * sinradangminusninety;
        double sy = poft * sinradangminusninety + doft * cosradangminusninety;
        return new Position(sx, sy);
    }

    public static double calcAmplification(Event event, double day) throws CalculationException {
        double amp;
        if (event.isSingleLens()) {
            amp = calcSingleLensAmplification(event, day);
        } else {
            amp = calcBinaryLensAmplification(event, day);
        }
        if (!Util.equals(event.getF(), 1.0)) {
            return blend(event.getF(), amp);
        }
        return amp;
    }

    public static double calcMagnitude(Event event, double day) throws CalculationException {
        return -2.5 * Math.log10(calcAmplification(event, day)) + event.getM0();
    }

    public static double blend(double f, double amplification) {
        return f * (amplification - 1.0) + 1.0;
    }

    /** Calculate amplification using Asada's method 
     * 
     * @param day   The absolute time of this data point
     * @param event The input event parameters
     * @return
     */
    public static double calcBinaryLensAmplification(Event event, double day) throws CalculationException {
        // first obtain source x- and y-position
        Position source = calcSourcePlanePosition(event, day);
        // Next get the Asada position parmeterization
        AsadaPositionParameters asadaPos = new AsadaPositionParameters(source, event);
        // now calculate the Asada image positions
        Position[] positions = calcImagePositionsUsingAsada(asadaPos, event.getQ());
        if (positions != null) {
            // calculate amplification from image positions
            double amplification = calcAmplificationFromImagePositions(event, positions);
            return amplification;
        } else {
            return Double.NaN;
        }
    }

    private static double calcLensEquationDeterminant(Event event, double x, double y) {
        double derx = 0.0;
        double dery = 0.0;
        // lens positions
        double[] px = new double[]{0.0, event.getA()};
        double[] py = new double[]{0.0, 0.0};
        double[] q = new double[]{1.0, event.getQ()};
        for (int i = 0; i < 2; i++) {
            double dx = px[i] - x;
            double dy = py[i] - y;
            double bottom = (dx) * (dx) + (dy) * (dy);
            bottom *= bottom;
            derx += (dx * dx - dy * dy) / bottom * q[i];
            dery += 2.0 * dx * dy * q[i] / bottom;
        }

        return 1.0 - (derx * derx + dery * dery);
    }

    private static double calcAmplificationFromImagePositions(Event event, Position[] positions) throws CalculationException {
        double totalamp = 0.0;
        for (int i = 0; i < positions.length; i++) {
            double temp = calcLensEquationDeterminant(event, positions[i].getX(), positions[i].getY());
            if (Util.equals(temp, 0.0)) {
                throw new CalculationException("Determinant zero");
            }
            totalamp += 1.0 / Math.abs(temp);
        }
        // a fix for missing images. These are probably due to bad Asada positions
//        if (positions.length < 3) {
//            totalamp += 3 - positions.length;
//        }
        return totalamp;
    }

    private static Position[] calcImagePositionsUsingAsada(AsadaPositionParameters asadaPos, double qForScaling) throws CalculationException {
        double r = asadaPos.getRho();
        double l = asadaPos.getL();
        double c = asadaPos.getC();
        double s = asadaPos.getS();
        double v = asadaPos.getNu();
        double[] a = new double[6];

        double r2 = r * r;
        double r3 = r2 * r;
        double r4 = r3 * r;

        double c2 = c * c;
        double c3 = c2 * c;
        double c4 = c3 * c;

        double s2 = s * s;
        double s3 = s2 * s;
        double s4 = s3 * s;

        double l2 = l * l;
        double l3 = l2 * l;
        double l4 = l3 * l;

        // start the coefficients proper
        a[5] = v * l * r3 * s3;
        a[4] = r2 * s2 + 2.0 * l * r3 * c * s2 - l2 * (2.0 * r2 * s2 - r4 * s2) - 2.0 * l3 * r3 * c * s2 + l4 * r2 * s2 - v * (5.0 * l * r3 * c * s2 - 4.0 * l2 * r2 * s2);
        a[3] = -2.0 * r2 * c * s - 4.0 * l * r3 * (c2 * s - s3) + l2 * (4.0 * r2 * c * s - 2.0 * r4 * c * s) + 4.0 * l3 * r3 * c2 * s - 2.0 * l4 * r2 * c * s + v * (l * (2.0 * r * s + 8.0 * r3 * c2 * s - 2.0 * r3 * s3) - 10.0 * l2 * r2 * c * s + 2.0 * l3 * r * s);
        a[2] = r2 + l * (2.0 * r3 * c3 - 10.0 * r3 * c * s2) + l2 * (-2.0 * r2 * c2 + 2.0 * r2 * s2 + r4) - 2.0 * l3 * r3 * c + l4 * r2 + v * (l * (-2.0 * r * c - 4.0 * r3 * c3 + 6.0 * r3 * c * s2) + 6.0 * l2 * r2 * c2 - 2.0 * l3 * r * c) + v * v * l2;
        a[1] = -2.0 * r2 * c * s + 8.0 * l * r3 * c2 * s - l2 * (4.0 * r2 * c * s + 2.0 * r4 * c * s) + 4.0 * l3 * r3 * c2 * s - 2.0 * l4 * r2 * c * s + v * (l * (2.0 * r * s - 4.0 * r3 * c2 * s + r3 * s3) - 2.0 * l2 * r2 * c * s + 2.0 * l3 * r * s);
        a[0] = r2 * c2 - 2.0 * l * r3 * c3 + l2 * (2.0 * r2 * c2 + r4 * c2) - 2.0 * l3 * r3 * c3 + l4 * r2 * c2 + v * (-l * (2.0 * r * c + r3 * c * s2) + 2.0 * l2 * r2 * c2 - 2.0 * l3 * r * c) + v * v * l2;

        // solve for roots
        int order = 5;
        double[] tanphi = new double[order];
        double[] rootsimag = new double[order];

        ComplexNumber[] roots = MathsEngineFactory.getMathsEngine(config.MATHS_ENGINE).polyroot(a);
        for (int i = 0; i < roots.length; ++i) {
            tanphi[i] = roots[i].real;
            rootsimag[i] = roots[i].imaginary;
        }

        // we should now have 5 values of tan(phi), one for each root!
        // If there are any imaginary roots, they are NOT REAL SOLUTIONS.
        // Hence, if inside a caustic, we'll have five real solutions, otherwise 3

        int realcount = 0;
        boolean[] realimage = new boolean[5];
        for (int i = 0; i < 5; ++i) {
            if (Util.equals(rootsimag[i], 0.0, 1.0e-4)) {
                realcount++;
                realimage[i] = true;
            } else {
                realimage[i] = false;
            }
        }
        if (realcount != 3 && realcount != 5) {
            return null;
        // by the way. LESS than 3 images are simply dealt with by adding 1 to amplification 
        // for each missing image
        }

        // now get r.cos(phi)
        double[] rcosphi = new double[5];
        for (int i = 0; i < 5; ++i) {
            if (realimage[i]) {
                double top = (l2 * r * c - l * r2 * c2 - v * l + r * c) * tanphi[i] * tanphi[i] - r * s * (l2 - 2 * l * r * c + 1.0) * tanphi[i] - l * r2 * s2;
                double bottom = r * (c * tanphi[i] - s) *
                        (r * s + 2 * (l - r * c) * tanphi[i] - r * s * tanphi[i] * tanphi[i]);
                rcosphi[i] = top / bottom;
            }
        }

        List<Position> imagePositions = new ArrayList<Position>();
        // don't forget to scale back to pvEvent scale!
        double scale = Math.sqrt(1.0 + qForScaling);
        for (int i = 0; i < 5; ++i) {
            if (realimage[i]) {
                imagePositions.add(new Position(rcosphi[i] * scale, rcosphi[i] * tanphi[i] * scale));
            }
        }
        return imagePositions.toArray(new Position[imagePositions.size()]);
    }
}
